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The stability of a quantized vortex state in Bose-Einstein condensation is examined within 
Bogoliubov theory for alkali atom gases confined in a harmonic potential under forced rotation. 
By solving the non-linear Bogoliubov equations coupled with the Gross-Pitavskii equation, the 
elementary excitations and the total energy of the systems are calculated as a function of rotation 
velocity. There are two distinct criteria of vortex stability; The position of the excitation energy 
levels relative to the condensate energy level yields the local stability criterion, and the total 
energy relative to that of the non-vortex state yields the global stability criterion. The vortex 
stability phase diagram in the rotation velocity vs the particle density of the system is obtained, 
allowing one to locate the appropriate region to observe the singly quantized vortex. 
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§1. Introduction 

There has been much attention focused on the Bose- 
Einstein condensation (BEC) experimentally and theo- 
retically|-j3iru:e it was realized in dilute alkali atom gases 
in 1995.BB0> The condensates have atom number typi- 
cally - O(I0 6 ) for 23 Na and 87 Rb, their BEC transition 
temperatures are in a range of (9(/iK), and theji .are usu- 
ally confined magnetically in a harmonic trap.Qu'EP 

In these systems, it is natural to expect a quan- 
tized vortex which has been already much studied in 
a supercoiiductor (Fermion) and 4 He (Boson) for a 
long time.lM While various microscopic theories for a 
vortex in Fermionic systems such as BCS-Gpr'kov or 
Bogoliubov-de Gennes theories are developed,© the cor- 
responding mean-field framework for a vortex in Bosonic 
systems origirialhL due to Bogoliubov has not fully ex- 
amined yet.oEj'E 1 ! Therefore, it is urgent to examine 
whether or not the widely used Bogoliubov theory, which 
has been quite successful so far for describing present 
BEC systems without a vortexJS 1 can also yield a stable 
vortex. 

In connection with the pres^L^djhrte^EC systems 
of alkali atom gases, theoreticaOBBBBBBB and 
experimental studies on a vortex have just started. So 
far there is no report to observe a vortex experimentally. 
One possible reason of the difficulty in producing a vor- 
tex is that gases are confined magnetically in a harmonic 
trap and this object is not so easy to manipulate it un- 
der rotation as compared with superflpd_lHe in a bucket. 
There are some theoretical proposalsOLLa) to realize the 
vortex state in this situation. The persistent rotation 
may be realized by shining moving laser on a large BEC 
system, which is examined theoretically by Marzlin et 

am 

Previously, we examined the stability problem of a vor- 
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tex in BEC systemsconfined by a rig id wal0 and by a 
harmonic potentialEU within the framework of Hartree- 
Fock Bogoliubov theory both at T=0 and T ^ 0. The 
conclusions we have reached are summarized as fol- 
lows: (1) At T=0 the vortex is unstable in the non- 
selfconsistent .Rogpliubov theory and the self-consistent 
Popov theory.ElEf (2) Above a certain temperature the 
vortex becomes stabilized by the presence of the in- 
creased non-condensate fraction localized at the vortex 
center, which effectively acts as. a pinning potential, pre- 
venting it from spiraling out.E3 (3) Even at T = 0, the 
external laser-induced potential at the_vortex center has 
a similar effect to stabilize the vortexO In these discus- 
sions, the stability of the vortex state ^.examined by the 
positive sign of the lowest eigenvalueO while its nega- 
tive sign means the eigenvalue becomes lower than the 
ground state energy. We call this stability "local sta- 
bility" . ThiL mipstinn of this type is still under lively 
discussions@BB@B@> 

As in 4 He svstem under persistent rotation around a 
symmetry axis,cP another stability criterion of the vortex 
state is to compare the total energy of the vortex system 
with that of the non- vortex system. We call this stability 
"global stability", considering the energy landscape of 
the configurational space. Here we study these two types 
of the vortex stability problem simultaneously, when the 
BEC system is kept in rotation at T = 0. 

It is expected that above a certain critical angular ve- 
locity figiobai, a vortex nucleates as in superfluid 4 He 
system or as in the lower critical field H c \ in a type II 
superconductor. As we will see in detail, the applied ro- 
tation tends to stabilize a vortex in general, therefore the 
above local stability criterion of the vortex formation is 
also satisfied even at T = 0, which is characterized by 
another angular velocity f2i OC ai- The purposes of this pa- 
per are to determine these characteristic angular veloc- 



ities; fiscal and figiobai and to examine the xe 
between them. Because the previous studyO' 
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corre- 
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sponds to Q — here, we will gain further understanding 
even for the local stability problem through this study 
for various fi's. 

To attack this problem, we solve the eigenvalue 
equation of the Bogoliubov theory, which is non- 
selfconsistent. If it gives a negative eigenvalue under 
a given condensate spatial profile determined by the 
associated no-linear Schrodinger type equation (Gross- 
Pitaevskii equation, or GP), the quantized vortex cannot 
be stable even for a self-consistent calculation. 

The arrangement of the paper is as follows: After giv- 
ing a brief description of the Bogoliubov theory in §2, 
the vortex stability problem in a system trapped har- 
monically is discussed in §3. The numerical procedure 
for solving the GP equation and the eigenvalue equation 
is the same as before, described in refs. ^l] and |2^. The 
final section is devoted to discussions and summary. 

§2. Formulation of Bogoliubov Equation Under 
Rotation 

2. 1 General formulation 

We start with a system of interacting Bosons. Two 
particles at r and r' interact with the potential g5(r — r') 
where g is a positive (repulsive) constant proportional to 
the s-wave scattering length a, namely g = 4:nh 2 a/m (m 
the particle mass). The system is kept under rotation 
with the angular velocity lo by an external force. The 
relevant hamiltonian H rot with the extra rotation term 
can be written as: 



H, 



H — lu ■ / r x p(r)dr 



(1) 



where 



H = / dr*t(r) <^ — - — + V(r) - fi \ *(r) 

+| [dr&(r)&(r)$(r)$(r), (2) 

p(r) is the momentum operator, /i is the chemical po- 
tential, and V(r) is the confining potential. 

Following the standard method, we decompose the 
field operator ^ as 



*(r) = -0(r) + <fr(r) 
0(r) = (#(r)). 



(3) 
(4) 



We substitute the above decomposition eq. (3) into eq. 
M) and ignore the higher order terms such as ip'ipip, 
tjr ^ ip , and ^ •ipxp terms. Then we rewrite the opera- 
tor ^(r) as 



(5) 



where r\ q (rj q ) is the annihilation (creation) operator and 
u q (r) and v q (r) are the wave functions, and the subscript 
q denotes the quantum number. 

2.2 Cylindrical system 

We consider a cylindrically symmetric system and a 
vortex line, if exists, passes through the center of a cylin- 
der, coinciding with the rotation axis. We use the cylin- 



drical coordinates: r = (r,8,z). The system is trapped 
radially by a harmonic potential 

1 



V(r) = -m(27rz/)V, 



(6) 



and is periodic along the z-axis whose length is L. We 
focus on the lowest eigenstates of the momentum along 
z-axis and the quantum number q in eq. (^) is written 
as (q r , qg), where q r = 1, 2, 3, • • ■ and qg = 0, ±1, ±2, • • •. 
The condensate wave function </>(r) is expressed as 



i,z) 



0(r)e i; 



(7) 



where <f>(r) is a real function and w is the winding num- 
ber. The case w = corresponds to a system without a 
vortex and w = l(w = 2) corresponds to a system with a 
singly (doubly) quantized vortex. The w > 2 case is not 
considered in this paper. The phases of u q (r) and v q (r) 
can be written as 



u,(r) =u q {r)e l{q <> +w 1 e 
v q (r) ^Vgiry^-^ 6 . 



(8) 
(9) 



The condition that the first order term in ip(r) of our 
hamiltonian vanish is 



2m [ dr 2 



1 d 

r dr 

1 2 



+V(r) + g<j> > (r) + Hlow <f>(r) = 



(10) 



This non-linear Schrodinger type equation is called the 
Gross-Pitaevskii equation. The condition that the hamil- 
tonian be diagonalized gives the coupled eigenvalue equa- 
tions for Uo(r), and v q (r) whose eigenvalue is e q : 

" I <r 1 u i n- (/(; ]- 1 

A* 



2m I dr 2 



H 

r dr 



+V(r) + 2g4> 2 {r) + huj(w + qg) u q {r) 

~gcj) 2 {r)v q {r) 
(w - qg) 2 ' 



e q u q {r) (11) 



2m 1 dr 2 r dr 

+V(r) + 2g<j) 2 (r) + huj{w - g e )J v q {r) 

-gcj) 2 {r)u q {r) = -e q v q {r\12) 

The normalization condition for 4>(r), u q (r), and v q (r) 
are 



/N 
\(j){r)\ 2 rdr = — , 



2tt / {u* p {r)u q {r) - v*(r)v q {r)} rdr = -6 P , 



(13) 
(14) 



where TV is the total particle number. 

The following symmetry relation should be noted: 

(u, v, e) at qg = n (v, u, —e) at qg = —n (15) 

We determine the signs of qg and e by the normaliza- 
tion condition eq. (^)- The set of the eqs. (10), (11), 
(12), (|l3|), and ( |l4| ) constitute Bogoliubov theory of our 
system. 
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When this Bogoliubov theory is extended to a finite 
temperature such as Popov approximation (see for exam- 
ple ref. |2l]), the eigenstates characterized by the eigen- 
functions u q and v q , and the eigenvalue e q are interpreted 
as the quasiparticles of the system which contribute to 
the total particle density as f(e)(\u q \ 2 + \v q \ 2 ) with f(e) 
being the Bose distribution function. 

The energy of the system is given by 



E 



(H rot ) + (iN 

£(m- 



E, 



2irL / \v q {r)\ 2 rdr. 



(16) 



Ea, = 2nL 



n ptjf d2 + 1 A _ ™1 

2m \ dr 2 r dr r 2 



-V(r)Wt 



<J 



rdr — wTiLoN. (17) 



Note that the coefficient of uo in eq. (17) is a constant. 
Although e's and u's are affected by the angular velocity 
u>, these contributions to E are negligibly small and we 
will use E,p as the total energy from now on. 

2. 3 Normalization 

The condensate has N particles in the length L and 
the area density is n z = N/L. The mass of a particle is 
m and the s-wave scattering length is a. The length r is 
measured as r' — (\j2mhv /K)r. We introduce the den- 
sity unit no = yj hv / g where g — A-Kfi 2 a/m. The density 
and energy are scaled by hq and hv respectively. There- 
fore, the normalized quantities are: 4>'{r') = -j=<j)(r), 

(if = rjzli- The angular velocity u> is also scaled by 2wv 
and the normalized rotation is Sfl ■ 



Equation (10) with eq. 
form 



1-KV " 

becomes in a dimensionless 



• d 2 


1 d 


w 2 ' 






dr' 2 ' 


r' dr' 






1 ,9 

r + 


0' 2 (/) 


+ wfl 


cf>'(r') = 0. 


(18) 



Equations (11) and (12) is transformed in a similar way. 
Equations (O) and (Oh are rewritten as 



'(r)\ 2 r'dr' = Aan z . (19) 



Wp (r')u'g(r') - v'*{r')v' q {r')} r'dr' = -^S p , q .{20) 



The energy E^ is rewritten as 



hvN Aan 
1 



-4>*'(r') 



V dr' 2 



1 d 

r' dr' 

r'dr' 



L 

2 



(21) 



Here 4an z and 4? are just dimensionless numbers. Since 
the parameter ^ does not change the following results 
within the present Bogoliubov framework, the system 
is characterized by the number an z and the normalized 
rotation f2 = 



§3. Two Kinds of Stability 

We solve the coupled equations; eqs. (10), (11), (12), 
and ( llil ) for the gas of 23 Na atoms trapped radially by 
a harmonic potential eq. ([|). The area density per unit 
length along the z-axis is chosen to be an z = 2.75 ~ 
137.5. (When the scattering length is a = 2.75nm, the 
area density n z varies from 0.1 to 5 (10 4 //im).) We use 
the normalized number f2 to indicate the angular velocity 
it). Q varies from to 1. 

Since the GP equation eq. (10) is detached from the 
rest of the equations, it is easy to know the spatial profile 
of the condensate cf> 2 (r) by solving the GP equation. Fig- 
ure 1 shows the typical density distribution of the system 
with and without the vortex along the radial direction. 
The doubly quantized vortex w — 2 case is also plotted 
for comparison. 
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Fig. 1. The particle number density distribution <f> 2 (r) along the 
radial direction. Dashed (dotted) line is the singly (doubly) 
quantized vortex case w = 1 (w = 2) while solid line corresponds 
to the non-vortex case. To draw this figure in actual units, we 
use the following parameters: the scattering length a = 2.75nm, 
m = 3.81 X 10 -26 kg, the radial trapping frequency v = 100Hz, 



and n. 



■■ 2 X 10 4 //im. 



3.1 Eigenvalues — local stability 

When the lowest eigenvalue of the eigenvalue equa- 
tions; eqs. (11) and (12) is negative in this formulation, 
the fundamental assumption of BEC is broken down as 
mentioned before. In this framework, the eigenstate with 
negative eigenvalue means that the system is unstable. 
Suppose that the negative eigenvalue appears in the vor- 
tex system while it dose not appear in the corresponding 
non-vortex system. Then the vortex state is unstable 
and-jprohibited against the non-vortex system. Dodd et 
alx3 first consider the lowest eigenstate problem in the 
vortex state under S7 = 0, suggesting the presence of the 
negative eigenvalue for a harmonic potential. 

We perform here the extensive computations for solv- 
ing the eigenvalue equations to know the lowest eigen- 
value under various external angular velocities Q. We 
show below the results of an z — 55 as an example. In 
Fig. 2, the lowest edge of the lowest eigenvalues, above 
which other eigenvalues are densely distributed, is plot- 
ted as a function of the quantum number qo for several 
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f2's. In the small ft cases the lowest eigenvalue occurs 
for the eigenstate characterized by the quantum num- 
ber qg = —1, while in the larger Q cases some positive 
eigenvalue of qg around qg ~ 10 becomes negative as ft 
increases. The eigenvalue at qg = does not change with 

n. 




-i 1 1 < 1 

-10 10 20 



Fig. 2. The lowest edge of the eigenvalues along qg for selected 
angular velocities f2's. Circles, filled circles, triangles, and filled 
triangles correspond to SJ = 0, 0.15, 0.3, and 0.45 respectively. 
The eigenvalues at qg = —1 move up while the eigenvalues in 
qg > move down as Q increases. 



We determine the critical Cl where the lowest e q 
changes its sign. The trace of the lowest eigenvalues as 
a function of VL is displayed in Fig. 3. The solid straight 
line is the eigenvalues at qg = — 1 and the dashed line is 
ones in qg > 0. The former becomes negative to positive 
at f^ ca i = 0.16085 while the latter becomes negative 
at fiscal = 0.44610. Therefore, the vortex state is lo- 
cally stabilized in the region 0.16085 < fi < 0.44610 for 
the present example of an z — 55. We call the stability 
defined by the sign of the lowest eigenvalue "local sta- 
bility". We define these critical fi's as 0} J ocal (= 0.16085) 
and OU cal (= 0.44610). 




Fig. 3. The trace of the lowest eigenvalues as a function of Q. 
Solid line corresponds to the eigenvalue of the (qg,q r ) = (— 1, 1) 
state and dashed line is other states with the lowest eigenvalue 
(q r = 1 and q e > 0). For the stable region 0.16085 < SI < 
0.44610 all the eigenvalues in the system are positive. 



3. 2 Interpretation of local stability 

The occurrence of this negative eigenvalue at qg = — 1 
may correspond to the 'spiraling out' of vortex predicted 
by RokhsarJia) Strictly speaking, the negative eigenval- 
ues only means that the system with a single vortex is 
unstable and prohibited. But the spatial variation of 
the eigenfunctions u q and v q suggests how the instabil- 
ity, which accompanies negative eigenvalue e qi actually 
occurs. 

When the sign of the energy of the lowest qg = — 1 
state changes from positive to negative at O = fij^ cal , 
this (excited) state must become the condensate state. 
Note that the original condensate was at E = 0. It is seen 
from Fig. 4 that this qg = — 1 state is localized around 
r = 0. Let us consider the transformation process of the 
condensate. We show the schematic diagram in Fig. 5. If 
the condensate is replaced by u(qg = —1) around r = 
whose effective angular momentum is w + qg = as seen 
from eq. (8) and still has the winding number w = 1 
at larger r region (Fig. 5), the vortex center as a phase 
singularity (A in Fig. 5) must exist at the border between 
the w + qg =0 region which is localized around r = 
and the w = 1 region which remains outside. The axial 
symmetry is broken and the vortex center slips out from 
r = 0. 

On the other hand, at il — nP Cf j, e of one particular 
state with qg > 0, say, qg = 10 becomes negative (see 
Fig. 2). The lower curves of Fig. 4 are the wavefunctions 
for the qg — 10 state. They are localized at the edge 
region of the condensate (j>. If the condensate has a large 
angular momentum w + qg ~ 10 at large r and has the 
angular momentum w = 1 around r — (see Fig. 6), this 
means that many phase singularities (i.e. vortices) must 
exist at the interface between the w + qg ~ 10 region at 
large r and the w = 1 region at small r. This transition 
stage will be followed by more complicated transition 
process. 
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Fig. 4. Wavefunctions u (solid line) and v (dashed line) of the 
lowest modes at qg = —1 (upper lines) and at qg = 10 (lower 
lines). The system have O = 0. These wavefunctions changes 
little even when f2 and the sign of e changes. The length unit is 
same as that of Fig. 1. 
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Fig. 5. The calculated wavefunctions <f>(r) and u(r) are showed 
in the right hand panel. They correspond to Figs. 1 and 4. The 
schematic figure of the phase change in the radial direction is 
drawn in the left hand panel. Because the new condensate with- 
out phase singularity appears at the center, the original vortex 
core (the triangle) must move outward. 




Fig. 6. The calculated wavefunctions <p(r) and u(r) are shown in 
the right hand panel. They correspond to Figs. 1 and 4. The 
schematic figure of the phase change in the radial direction is 
drawn in the left hand panel. As the new condensate at large r 
has large (~ 10 X 2tt here) phase change, many vortex cores (the 
triangle) must appear. 



3.3 Total energies — global stability 

The comparison of the two energies in the systems with 
and without a vortex is another measure of the vortex 
stability. The comparison indicates the global stability 
of the vortex state relative to the non-vortex state. We 
calculate the critical 51 in which the energy of the system 
with the quantized vortex is equal to the system without 
a vortex under the same condition. We ignore the energy 
coming from the non-condensate part because the energy 
is 100 times smaller than the energy of the condensate 
in the present situation at T = 0. 

In the non-vortex state, the energy E w —o is not af- 
fected by the rotation 51, as is seen from eq. (17). The 
energy of the system without the vortex E w —o, the en- 
ergy of the singly quantized vortex system E w —i (51), and 
the energy of system with the doubly quantized vortex 
E w= 2(£l) are evaluated from eq. (21) as 

E W=1 (Q) - E W=Q = hvN {0.22215 - Q) (22) 



E w=2 (il) - E w=1 {fl) = huN {0.42071 - 51). (23) 

Therefore, the critical 51 g i b a i where the two energies are 
equal; E w= q — E w= i(iY) and E w= i(fl) = E W —2(Q) are 
given by ^ lobal = 0.22215 and ^ lohal = 0.42071 for 
an z — 55. We call the stability defined by the energy 
crossing "global stability" . 

3-4 Two stabilities 

While the vortex stability in 4 He systems is usually 
discussed in terms of the total energy (global stability 
in this paper), the .stability in the present dilute BEC 
systems is discussecO using the lowest eigenvalue (local 
stability in this paper). These two stabilities are deter- 
mined here and we can compare these quantitatively for 
a wide range of the parameter values of an z . 

The local stability is satisfied for 0.16085 < 51 < 
0.44610 and the global stability is given by 0.22215 < 
51 < 0.42071 when an z = 55. In the substantially wide 
51 region the vortex state of the rotating BEC systems is 
stable locally and globally. 

As mentioned previously, the system is characterized 
by the number an z and the rotation 51. In Fig. 7 we 
show four critical 51's as a function of an z . It is seen 
from Fig. 7 that as 51 increases, (1) the singly quantized 
vortex becomes stable locally first, (2) the energy of this 
vortex state is lower than that in the corresponding non- 
vortex state and the singly quantized vortex nucleates 
in the system, corresponding to (B) in Fig. 7, and (3) 
this single vortex state is kept stable up to either 51g lobal 
or 5l[ ocal . These divide the parameter space into three 
regions: the vortex unstable region (A), the single vortex 
stable region (B) and the single vortex unstable region 
(C) where many vortices may appear. 





(C) 
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~~cfi " — ~ — 

^ ^global 


11 ^local 


(A) 





50 100 140 

an z 

Fig. 7. Phase diagram in the Q vs an z plane, showing the three 
characteristic regions: the vortex unstable region (A), the single 
vortex stable region (B) and the single vortex unstable region 
(C) where many vortices may appear. Solid lines show n g i b a i 
and dotted lines show f2i oca j. Upper two lines and lower two 
lines are f2 u and f2 L respectively. 



§4. Conclusion and Discussions 

In order to obtain more insights into the vortex stabil- 
ity problem in the Bose-Einstein condensation of alkali 
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atom gases confined in a harmonic potential, we have 
extended our previous workEUEif to the case where the 
cylindrical system is under forced rotation. 

Our calculation is done within the framework of the 
non-selfconsistent Bogoliubov theory and yields the vor- 
tex stability phase diagram shown in Fig. 7. The stabil- 
ity of the vortex state is examined by the two different 
ways, that is, the localj-stability and the global stabil- 
ity. The previous cases jEjEj which are consistent with 
the present results, correspond to the O = axis in this 
figure where the vortex is intrinsically unstable. This in- 
stability is seen to occupy a finite region, rather than the 
isolated line confined at — for various densities n z 
and the scattering lengths a. 

We have examined not only the above intrinsic local 
stability, but also the global stability of the single vor- 
tex relative to the non-vortex state. This allows us to 
draw the whole perspective for the vortex stability prob- 
lem. As is seen from Fig. 7, this global stability bor- 
derline Og lobal corresponding to the vortex nucleation is 
always situated above the intrinsic stability borderline 
^focai- This means that for any an z values investigated 
the singly quantized vortex, which is nucleated, can exist 
as a (locally) stable object under rotation at T=0. It is to 
be noted that the present results in the non-selfconsistent 
calculations are applicable even for selfconsistent calcu- 
lations such as the Popov approximation. 

We have performed a similar study for the rigid wall 
case instead of the present harmonic potential case. Our 
data show the same vortex stability phase diagram for 
the system sizes up to 80 x /coherence length), confirm- 
ing the previous conclusionEr that the intrinsic vortex 
unstable region exists at a lower £1 region. 

Together with the previous finite temperature calcu- 
lations,^ the present calculation concludes that alkali 
atom Bose gases in BEC can sustain and exhibit the sta- 
ble vortex in appropriate temperatures and appropriate 
rotations. 
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